The shower approach in the simulation of ion scattering from solids 



V.A. KhodyrevQ R. Andrzejewski, A. Rivera, D.O. Boerma, and J.E. PrietcQ 

Centro de Microanalysis de Materiales and Instituto Universitario "Nicolas Cabrera", 
Universidad Autonoma de Madrid, E-28049 Madrid, Spain 

For basic studies of ion-solid interactions as well as for the retrieving of reliable structural in- 
formation from ion scattering experiments, accurate simulations of swift ion-solid interactions are 
essential. This, however, faces a fundamental problem: due to the small cross sections for scat- 
tering, the probability of an incoming particle to reach the detector is always small. As a result, 
programs based on the "brute force" approach where ion trajectories are calculated taking into ac- 
count the large diversity of histories of collisions, require extremely long computation times to reach 
acceptable statistics. In this paper, a new simulation approach is proposed which combines an exact 
treatment of the interactions within the binary collision model with a high computational efficiency. 
The main idea is to use the strategy of importance sampling in the simulation of atomic thermal 
displacements: the positions resulting in scattering events close to the direction of the detector are 
sampled with increased probability. As a result, the detector is illuminated by intensive showers of 
scattered ions, even if multiple scattering is important in the way out of the sample volume. The 
developed computer code provides the possibility to address wide classes of simulation problems. 
It is able to treat complex crystal and amorphous structures, to use different models of ion-atom 
interaction (potential, inelastic energy loss, multiple scattering on electrons). The program provides 
information such as depth profiles of close collision yield, energy spectra and 2D angular distribu- 
tions of scattered and recoiled ions. To illustrate these features we present the results of simulations 
for several, widely different conditions and compare them with experimental results. It is argued 
that, in fact for the first time, this method provides a possibility to reliably treat the rare events of 
plural scattering, the most intricate problem for the theoretical description as well as for computer 
simulations of ion-solid interactions. 

PACS numbers: 02.70.-c, 61.85.-pp, 34.35.+a 



I. INTRODUCTION 



The classical description in the binary collision approx- 
imation is well applicable for the description of ion-solid 
interactions in the energy range above ~1 keV; the mul- 
tiple interaction with two or more atoms, if significant, 
may be treated in a simple first-order approximation. In 
the computer simulation of phenomena related to ion- 
solid interactions, this provides the possibility of an effi- 
cient reproduction of many experimental conditions play- 
ing with the parameters of the interaction model. Due 
to the importance of this subject great efforts have been 
performed towards the development of efficient simula- 
tion programs [l[ . To illustrate the central role that com- 
puter simulations have played in the field, it will suffice 
to mention that one of the most prominent phenomena, 
the channeling of ions in crystals, was first observed in 
simulation results Q. In this context, one should refer 
also to the paper of Barrett [i| which provides an im- 
portant supplement to the theory of ion-crystal interac- 
tion allowing to address some aspects of the phenomenon 
which arc difficult to treat theoretically. Currently, with 
ion beams being widely used as a precision tool for the 
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analysis and modification of materials, there has been 
an increase of the level of demands to simulation algo- 
rithms. One example is the use of low-energy ion beams 
for surface structure determinations. The classical pic- 
ture of scattering together with the effects of blocking 
of scattering on atomic pairs form the basis for obtain- 
ing detailed information about the atom locations. The 
only way to extract this information is the comparison of 
intensities measured for different scattering geometries 
with the results of simulations for many trial structures. 
However, a fundamental problem occurs here due to the 
insufficient efficiency of existing algorithms. Due to the 
small scattering cross sections, the direct simulation by 
calculation of individual ion trajectories (the program 
MARLOWE [3] is the most developed code of this type) 
is often not practicable. To understand this one has 
to keep in mind that the experimental procedure typi- 
cally consists in the measurement of angular scans with 
a small-aperture energy detector. To acquire the nec- 
essary statistics in a measured spectrum, the required 
number of ions in the beam amounts to ~ 10 9 for low- 
energy (keV range) scattering of heavy ions like Ne or 
Ar and to ~ 10 13 for scattering of ions with energies in 
the MeV range. Specially in the latter case it is out of 
the question to calculate the histories of motion of such a 
number of particles directly. For this reason, the conclu- 
sion is reached [l[ that even the power of supercomputers 
is by far not enough to perform such simulations. 

An overview of the existing approaches to this prob- 
lem shows that two main ideas are used depending on 
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the ion energy. In the simulation of backscattering of 
high-energy channeled ions, one can rely on the single- 
scattering model assuming that the motion of the ion 
in the outgoing path can be described by a straight-line 
trajectory. This allows to avoid a precise description of 
this segment of the ion trajectory. As a result, each 
ion contributes to the statistics according to the prob- 
ability of close collisions along the ingoing path. The 
inverse (blocking) condition can be treated analogously 
(see Ref. [f|, for example). This model is well tractable 
and serves as the basis for a large number of algorithms 
which have been developed (see Refs. [1, @, S| for the 
most widely used). 

This algorithm does not work, however, if one is inter- 
ested to simulate the channeling-blocking conditions or 
the scattering of low-energy and/or heavy ions. In these 
cases, significant influence of crystal structure and/or 
strong multiple and plural scattering result in the ne- 
cessity to reproduce the whole history of motion of in- 
dividual ions inside the sample. As a result, the single- 
scattering algorithm needs a serious modification: now, 
to avoid the calculation of outgoing trajectories, we 
need to know the probability for an ion to end up in 
the detector for every position of the atom of the con- 
sidered lattice site. A solution with emphasis for ap- 
plications in the medium-energy ion scattering (MEIS) 
case (~ 100 keV light ions) was proposed in Ref. [9| 
(the computer program VEGAS). According to the time- 
reversibility rule [lOj, the mentioned probability can be 
represented by the flux of ions at the same point in the 
case that the beam enters, inversely, from the detector 
side. As a result, the simulation consists in the calcula- 
tion of two fluxes, the direct flux and the inverse flux, 
followed by their convolution with the probability that 
the scattering atom occupies the required position. In 
this approach the electronic and nuclear energy losses 
are neglected, thus its applicability is guaranteed only at 
small depths of ion penetration and for relatively light 
ions. 

Due to purely technical difficulties (the main one is 
the need to store large sets of data for the fluxes), this 
method for medium energies is always used in a simpli- 
fied form. Namely, the energy and angular variables are 
not considered during the convolution of fluxes. It is as- 
sumed, in fact, that there is a one-to-one correspondence 
between energy and depth and the scattering angle of the 
collision is taken to be the same as the angle between the 
direction of the incoming beam and the direction to the 
detector. Thus, such approach considers mainly the ef- 
fect of the spatial distribution of fluxes and cannot claim 
to achieve a detailed description of energy spectra and 
multiple (plural) scattering effects. 

In the case of low-energy (E ~ 1 keV) ion scatter- 
ing (LEIS) where only a small number of atomic layers 
contribute to the detected yield, the full version of this 
"reversing approach" seems to be possible ll| . The only 
additional problem arising is that with the higher dimen- 
sion of phase space the accumulated data arrays become 



very dilute. As a result, the convolution of the fluxes re- 
quires some special procedure. For this purpose, it was 
proposed to sum pairs of crossing trajectories of incoming 
and (time-reversed) outgoing fluxes with weights deter- 
mined by the scattering cross section. It is clear that 
some tolerances in coordinates and energy have to be as- 
sumed in order to identify crossings among the finite set 
of calculated trajectories; the associated inaccuracies are 
equivalent to those present in the case of the numerical 
convolution of both fluxes. One might expect, however, 
that allowing for an admissible level of accuracy, the ef- 
ficiency of the simulation can be considerably increased 
as compared to the direct approach. Indeed, a satisfac- 
tory description of experimental results for several LEIS 
conditions was found as a result of the application of this 
approach fl2| . 

In the present paper a new simulation approach is pro- 
posed. The developed computer code has the name TRIC 
(Transport of Ions in Crystals). In this method, complete 
histories of ion collisions inside the sample are calculated 
as in the direct simulation approaches. To increase the 
efficiency we use a specific advance of the Monte Carlo 
method, the strategy of "importance sampling" [l3[: the 
distribution of thermal displacements of crystal atoms is 
sampled giving preference to those positions which result 
in scattering in a direction close to that of the detec- 
tor. The actual probabilities of such displacements are 
accounted for in the accumulated statistics by the use 
of appropriate weighing factors. This procedure results 
in an illumination of the detector by intensive "showers" 
of scattering events. Additionally, in order to achieve a 
smoother statistics of rare events due to plural scatter- 
ing, the strategy of "stratified sampling" [14[ is applied: 
the main fraction of close collisions which results in scat- 
tering to different directions are also sampled more often 
and these trajectories are also allowed to send showers 
towards the detector. The same strategy is used for a 
preferable choice of the initial coordinates of the ions en- 
tering into the crystal among those resulting in higher 
efficiency in the accumulated statistics. 

The paper is organized as follows. We present in 
Sect. II a detailed description of the simulation proce- 
dure. In Sect. Ill the efficiency and quality of the proce- 
dure are illustrated by the description of several exper- 
imental results under different conditions, selected also 
based on general interest. Since the reversing approach, 
closely related with our proposed method, is in fact the 
only alternative claiming a similar efficiency we pay spe- 
cial attention in Sect. IV to a careful analysis of its basis; 
the main intention there is to recognize clearly the rela- 
tive advantages of the two approaches. Some concluding 
remarks are made in Sect. V. 



II. THE SIMULATION PROCEDURE 

The main idea that allows the increase of efficiency of 
the simulation is based on the fact that, due to ther- 
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FIG. 1: (color online) The principle of selection of the "hot" 
region of atomic displacements used in the shower generation. 
The density of the Gaussian distribution of atomic displace- 
ments is shown by the contrasted gray while the "hot" region 
is indicated as the tube enclosing the region of impact pa- 
rameters for scatterings into the chosen angular cone (left). 
Analogously, the picture at the right panel illustrates the case 
of showers of recoiled atoms 

mal motion, each crystal atom can be found in different 
positions when the incoming ion crosses the considered 
lattice site. In the case that the atom is located close 
to the ion trajectory, a significant change in the ion mo- 
tion can take place. Events of scattering into directions 
close to the direction to the detector are obviously of 
special interest and it appears desirable to treat them 
selectively. This procedure should increase the number 
of relevant events at least in cases where scattering does 
not strongly perturb the further motion of the ion. 

To implement this simple idea we use the following 
procedure. First, the impact parameter bi which corre- 
sponds to scattering from the current direction of mo- 
tion into the direction to the center of the detector (see 
the left-hand side of Fig. [T]) is determined (i denotes the 
number of the current lattice site). Knowing this value, 
we can define the line (in the scattering plane) of atom 
positions necessary for such scattering. Then, assuming 
that the relevant scattering directions lie within a cone 
of width AO, we can define the corresponding region of 
impact parameters whose half-widths in the scattering 
plane, A6||, and in the perpendicular direction, Ab±, can 
be estimated as 

A6|| « |cZe/d6| -1 A6, Ab ± w (6/ sin 6) AG. (1) 

In this way also a fraction of the distribution profile of 
atom displacements is separated, the "hot" region. Then, 
one or several (n) atom positions can be drawn from this 
region treating it as a part of the total (Gaussian) dis- 
tribution of atomic displacements; as a result a shower 
of ion trajectories is directed towards the detector. The 
remaining part of the distribution is used to continue the 
trajectory of the "primary" ion. In such a way all possi- 
ble displacements of the atom are sampled. Notice that 
the width of the shower AG is assumed to be significantly 
larger then the width of the detector 50. 



In the accumulated statistics, events due to ions in the 
showers ( "secondary" ions) must contribute with weights 
Wi given by the probability for the atom to be displaced 
into the "hot" region. To calculate them we notice first 
that a probability should be adscribed also to the pri- 
mary ion itself (Wi before the i-th collision) which de- 
creases with depth as a result of emission of the secondary 
ions, (trajectory "degradation"). Then the weights Wi 
and Wi+i after the collision are updated as 

Wl = PiWi/n, W i+ i = (1 - PJWu 

where Pj is the integral fraction of the probability of 
the "hot" region. The sum of probabilities is conserved, 
Wi+i + nwi = Wi, as necessary. It is worthwhile to 
note also that this approach provides absolute values of 
the scattering yield: the expectation value for a certain 
energy-angular range is equal to that obtained in a di- 
rect simulation with the same number of ions sent to the 
sample. 

One can easily recognize the similarity of this approach 
with the strategy of "importance sampling" used in the 
Monte Carlo numerical integration [13[ where the val- 
ues of the integration variable are sampled more densely 
in the regions that give the highest contribution to the 
integral. This similarity is not accidental because the 
scattering yield, as an average over the ion initial con- 
ditions and the thermal displacement of atoms, is, in 
fact, determined by an integral over a multidimensional 
space. Due to the inherent complexity of this integral, 
the importance-sampling strategy can be used only on 
an intuitive basis as described above. It is clear that this 
strategy should be effective at high energies where the 
single-scattering model describes well the process. At 
low energies the shower approach can also be efficient. 
In this case the showers emitted at large depths arc sig- 
nificantly attenuated due to strong re-scattering on the 
way out of the crystal. This is a disadvantage but it is 
largely compensated by the fact that the primary ions 
themselves can leave with large probability the sample 
volume and the numerous accompanying showers from 
the top layers produce an intensive flux in the direction 
of the detector. It can be said that, in this case, the 
main part of the procedure is the direct simulation of ion 
trajectories and the showers are used only as an improve- 
ment at the last stage. These arguments show that the 
most serious difficulties in the application of the shower 
approach can be met in the simulation for intermediate 
ion energies where the rare plural-scattering events result 
in violent fluctuations (see examples in the next section). 

As shown at the right-hand side of Fig. [TJ a similar 
procedure can be used for the simulation of the yield of 
recoiled atoms. When an ion crosses a lattice site occu- 
pied by an atom of the considered species, a shower of 
recoils is emitted. The "hot" region selects now the atom 
displacements which result in emission of recoiled atoms 
within the chosen cone. One additional recoil atom is sent 
by sampling the remaining part of the Gaussian distribu- 
tion of atomic positions. This recoil, in turn, is allowed 



4 



to produce new recoils, as in the case of the primary ion, 
and also scattering showers. The trajectory of the im- 
pinging primary ion is followed further with the atom 
displacement drawn according to the total Gaussian dis- 
tribution. To describe the whole cascade we consider also 
recoils produced by the ions in showers (not allowed to 
produce new showers). The result of these many possi- 
bilities is a strongly developed tree of the cascade. Since 
all particles in the cascades have similar histories of col- 
lisions, the treatment using a recursive algorithm turns 
out to be very efficient. 

Violent fluctuations in the accumulation of scattering 
events can have different causes. A first important point 
in this respect is that, when the showers are generated 
at every lattice site crossed by the primary ion, this ion 
itself cannot scatter into the detector direction. This ex- 
cludes the possibility of detection of events with weights 
Wi which are ordinarily much larger than the weights of 
the secondary ions Wi. Thus, the displayed fluctuations 
are entirely due to the dispersion of the values of Wi . 

If the number of ions in one shower n is fixed, the value 
of the weight Wi for a certain scattering angle O depends 
on the position of the "hot" region with respect to the 
center of the Gaussian distribution of thermal displace- 
ments. Since the width of this distribution is rather small 
compared to interatomic distances, variations of several 
orders of magnitude in the values of Wi are possible. To 
improve upon this situation, it is reasonable, instead of 
fixing n, to fix the value of the weight of one event wq, 
treating accordingly the number of ions in a shower as 
variable, rii = Wi/wo (the fractional part is treated as 
the probability for sending one additional ion). When 
Wi 3> wo the resulting effect is expressed as an increase 
of the number of detected events with a smaller weight 
wo from intense showers (distributed to some extent over 
angles and energies) instead of one event with a large 
weight. An advantage in the inverse case, Wi <C wq, 
is that the load on the computer due to the simulation 
of non-significant events is avoided. In addition, the 
discrete counts of such "quanta of probability" closely 
mimics the aspect of experimental data. For orientation, 
a reasonable value of the quanta wo can be estimated 
defining n in the input as the number of ions in the most 
intense shower (wi = max). 

In the case of crystalline structure of the sample one 
has a further possibility to increase the number of de- 
tected events. Clearly, trajectories of ions entering into 
the crystal at different points of the surface can produce 
significantly different contributions to the useful statis- 
tics. Thus, in the sampling of the initial coordinates it 
seems reasonable to choose points with higher efficiency 
more often. In the program TRIC this is organized as 
a self-learning procedure: the distribution of initial co- 
ordinates, taken initially as uniform over the unit cell of 
the crystal surface, is deformed in the course of the sim- 
ulation according to the accumulated information on the 
efficiency. As necessary, the initial value of the weight 
Wo is taken to be inversely proportional to the density 



of the simulated distribution. Effectively, such modifi- 
cation corresponds to a larger number of impinging ions 
and this results in faster (up to an order of magnitude) 
accumulation of statistics. 

Another source of fluctuations is the plural scattering. 
In order to develop a picture of how these fluctuations 
arise we should mention the following. An appropriate 
choice of the width of the shower AO is determined by 
the demand that the shower should be wide enough to 
encompass the whole profile of multiple scattering in the 
outgoing path. One may believe in this case that the 
transport inwards and the transport outwards the central 
region (of width SO) are properly balanced in the shower. 
However, at low energies and/or for heavy ions, this con- 
dition is difficult to fulfill because the shower should be 
taken very wide and, as a result, only a small fraction 
of the ions in the shower reaches our small- aperture de- 
tector. But, if the above condition is not fulfilled, the 
reduction of ion flux in the shower cone due to diffusion 
outward is to be compensated by the plurally scattered 
primary ions. As this takes place, the main contribu- 
tion comes from primary ions scattered accidentally to 
directions close to the cone mantle. Since the probability 
of further scattering into the cone is large for these ions 
(small \dQ/db\ in (JTJ) ) they produce rare but very intense 
showers of secondary ions with almost equal energies. As 
a result the accumulated energy spectrum is disturbed 
by splashes in some channels. 

We solve this problem in a similar way as described 
above: the plural scattering events are sampled more 
frequently than they happen in reality. For this goal a 
second cone is considered, coaxial with the first and sig- 
nificantly wider than it. Wide showers are produced by 
the primary ion inside this cone in addition to the show- 
ers in the internal cone. Then, ions of the outer cone 
are allowed to send showers into the internal cone, simi- 
larly to the primary ion. As a result, our goal is achieved 
because the weights of ions in the outer cone are small 
compared with the weight of the primary ion Wi and the 
number of showers sent to the internal cone is now large 
to smooth sufficiently the accumulated energy spectrum. 
In principle, a whole hierarchy of such nested cones or 
even some smooth deformation of the density of sampling 
of atom displacements could be organized. It turns out, 
however, that in many cases the implementation of the 
above two-step approach solves practically all problems. 
This is illustrated in the next section (Fig. [7^). 

To estimate the efficiency of the shower approach com- 
pared to the direct simulation we should consider first 
the values of the weights of the accumulated events wq ■ 
Since the weight of an event in the direct simulation is 
unity, the number of impinging ions in the latter should 
be larger approximately by a factor 1/wq to achieve the 
same level of fluctuations in the statistics. Of course, 
this cannot serve as an accurate estimate of relative ef- 
ficiency, the calculation of the history of motion with 
shower emission is more extensive while, on the other 
hand, the calculated exiting flux is concentrated near 
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the relevant region. Nevertheless, the direct compari- 
son of two approaches, where possible, shows that the 
mentioned factor reflects well the situation. The value of 
this factor depends on ion energy and is typically ~ 10~ 3 
at keV energies while at hundreds of keV it is ~ 10~ 7 . 
Therefore, the advantage of the shower approach in this 
respect is enormous. A typical time for simulation of one 
energy spectrum using an ordinary PC is several minutes, 
thus the proposed approach is capable to treat problems 
practically not solvable by the direct approach. Note that 
both approaches are based on the same, rather general, 
model of ion-solid interaction. 

It is worthwhile to note also that the sending of show- 
ers aiming to the detector mainly results in misses. This 
unavoidable drawback, when one is interesting only in 
the yield of a small-aperture detector, turns out to be 
an advantage if one needs to simulate the 2D angular 
distribution of the scattered ions. Such possibility is il- 
lustrated in the following section. In fact, we are capable 
to calculate with the present approach 3D-distributions 
including the energy scale. 

The implementation of this algorithm as a FORTRAN 
program includes all important elements of the binary 
collision model [l[ . A wide possibility to choose the type 
of the interaction potential is provided. The scattering 
integrals for binary collisions as functions of impact pa- 
rameter and energy are tabulated at a preliminary stage 
of the calculation and used afterwards by interpolating 
with splines. Additionally, the impact-parameter depen- 
dence of inelastic energy losses is considered. To account 
for the simultaneous interaction with two or more atoms 
we sum the deviations in ion motion due to the interac- 
tion with individual atoms. In principle, all this provides 
the possibility to perform simulations for energies in a 
wide range. At high energies the procedure of generation 
of showers could encounter difficulties due to restrictions 
in the accuracy of computer calculations. The reason 
is that the width of the "hot" region in the transverse 
plane becomes so small that an accurate transformation 
of impact parameter to scattering angle may be practi- 
cally impossible. For this reason, if this width becomes 
smaller than 0.1 of the thermal vibration amplitude u\, 
the procedure of shower generation is inverted: the scat- 
tering angle is sampled within the shower cone and then, 
in reverse order, the impact parameter is calculated. As 
a result, our approach can claim to represent an improve- 
ment of the single-scattering model as it additionally ac- 
counts for the multiple scattering on the outgoing path. 

Materials of any crystal structure, including com- 
pounds, with rectangular unit cells can be described in 
the input. On a lower level, however, the program works 
with a description in terms of the Wigner-Setz cell and, 
in principle, this provides the possibility to treat also 
wider classes of structures. An amorphous structure can 
be simulated by a random rotation of the crystal lattice 
after each collision. The surface of the sample is de- 
fined as an imaginary plane appropriately located with 
respect to the crystal lattice. If the detector position is 



defined to be at the back side of the sample, a simulation 
of transmission through a crystal slab of a given thick- 
ness is performed. As output data, the depth profiles 
of close collisions, the energy spectra and the 2D angu- 
lar distributions of scattered ions or recoiled atoms of a 
certain species are calculated. Finally, the procedure is 
automated to calculate angular scans with step-wise ro- 
tation of the target, the detector or the beam direction 
around a certain axis. The FORTRAN code of the pro- 
gram TRIC is supplied with a graphical user interface 
allowing to comfortably supply the input data, to run 
the calculations and to inspect the simulation results. 



III. PROGRAM TESTS AND APPLICATIONS 

In this section we present examples of applications of 
our program based on the shower approach to simulations 
under different experimental conditions. The first one 
concerns the scattering of 5 keV Ar + ions from the (100) 
surface of a Fe4N crystal as studied experimentally in 
Ref. [lj| . The structure of this crystal can be considered 
as fee Fe with an additional N atom located at the center 
of the unit cell. The topmost atomic layer at the surface 
contains nitrogen and, depending on the growth condi- 
tions, can exist in two types of surface reconstructions, c4 
and 4pg. Figure shows simulated angular scans for the 
c4 surface differing from the bulk-terminated one only by 
an outward displacement (0.23 A) of the nitrogen atoms. 
In these scans, the crystal is rotated around the surface 
normal, which is coplanar with the beam and the detec- 
tor directed respectively at 42° and 12° from the surface 
at opposite sides of the normal. The yield of scattered 
Ar as well as of Fe and N recoils were calculated. In 
Fig. [3] the energy spectrum of scattered Ar is shown for 
the orientation indicated in Fig. [2] by an arrow. Also 
shown is the angular distribution of scattered ions. The 
latter illustrates the blocking of scattering from the top- 
layer Fe atoms by re-scattering from nearby atoms; this 
largely explains the strong variation of the yield in the 
angular scan. These variations are strongest when ei- 
ther the beam or the detector are positioned close to 
the scattering "horizon" which is formed by the reflec- 
tion from the surface as a whole (very small intensity 
at the bottom part of the angular distribution). In the 
case shown in Fig. [3J the detector is at the border of the 
shadow cone. The shape of the energy spectrum as an 
isolated peak is due to the blocking of the particles scat- 
tered from atoms of deep atomic layers by atoms of the 
top layer. This shape of the spectra favors the choice of 
the energy window for the performance of angular scans. 
The double-humped structure is explained by the effec- 
tive competition of sin gle- and double-scattering from Fe 
atoms of the top layer [16| . 

As seen in Fig. [2 the experimental results are well re- 
produced in the simulation both for scattered ions and 
for recoils of the two atomic species. Note that, in or- 
der to achieve this agreement, the screening radii in the 
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FIG. 2: Comparison of simulated angular scans for the scat- 
tering of 5 keV Ar + ions from the surface of a Fe4N(100) crys- 
tal with the experimental results from Ref.pjl (top panel). 
The yield of recoiled Fe and N atoms is shown in the center 
and lower panels, respectively. See text for more details. 



used "universal" potentials [TlJ were properly reduced, 
as commonly done for the description of scattering of 
heavy ions at low energies. With the same potential the 
data for the more complex reconstruction 4pg are also 
described with the same quality. In general, the present 
results demonstrate that the quality of the description of 
the experimental data is similar to that provided by the 
code MATCH [l5|; in fact, the results of both simulation 
can hardly be distinguished when plotted together. 

The second example, shown in Fig. 0^,, corresponds 
also to a LEIS experiment, now of 3 keV Ne + scattering 
from a polycrystalline Cu sample. Under these condi- 
tions an intriguing peak is observed [l8| in the experi- 
mental spectrum at the energy corresponding to a single 
Ne-Cu collision. Here the beam incidence was along the 
surface normal and the scattering angle was 129°. For 
these strong-scattering conditions the code MARLOWE 
is capable to reproduce the general shape of the spec- 
trum including the surface peak. Our simulations shown 
in the same picture also reproduce the spectrum shape; 
some differences with the MARLOWE simulation at low 
energies may be due to a difference in the used poten- 
tials of ion-atom interaction or in the treatment of in- 
clastic energy losses. In general, the shape of the spec- 
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FIG. 3: Energy spectrum of Ar + ions scattered from the (100) 
surface of a Fe4N crystal for the orientation indicated in Fig. [2] 
by an arrow. The lower part of the figure shows the angular 
distribution of scattered ions. Results are shown using a gray 
scale with darker representing higher intensity. 



trum under these conditions differs considerably from the 
spectrum predicted by the single-scattering model (also 
shown in the figure) . The authors of Ref. [18[ associated 
the appearance of this peak with the onset of plural and 
multiple scattering in the deeper layers. First, due to 
the the large variations of kinematic energy losses, the 
plural collisions in scattering from deeper layers yield a 
broad scattering spectrum (in particular, this explains 
the presence of the high-energy tail) while (single) scat- 
tering events from the top layers form the narrow peak. 
On the other hand, the latter is also favored by the rela- 
tively small inelastic energy losses: independently of the 
depth of scattering, all single-scattering events have simi- 
lar energies at the exit. To demonstrate the decisive role 
of strong localization of energy losses at small impact 
parameters, we present in Fig. 0b the results of a simula- 
tion where both nuclear and electronic energy losses are 
replaced by equivalent continuous stopping forces (thin 
solid curve). As seen, with such description of energy 
losses, the spectrum changes drastically. First, the sur- 
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FIG. 4: (a) Comparison of our simulation results for 3 keV 
Ne + ions scattering from amorphous copper with experimen- 
tal data [l8|| and results of simulation with the MARLOWE 
code. The simulation results were folded with a gaussian dis- 
tribution in order to simulate the experimental energy resolu- 
tion of 40 eV. The dashed curve corresponds to a calculation 
in the single-scattering approximation. Panel (b) shows sim- 
ulations performed by applying additional approximations. 
The dashed curve corresponds to a simulation using the single 
scattering approximation with impact-parameter-dependent 
energy losses while the thin solid curve was calculated con- 
sidering the energy losses as the result of a uniform stopping 
force. 



face peak disappears since single scattering events from 
different depths have now different energies. And second, 
the variation of kinematic energy loss in plural scattering 
is now not effective and, as a result, the shape of the spec- 
trum on both sides of the surface peak is also strongly 
modified. On the other hand, in the shower approach we 
can check also the role of plural and multiple scattering. 
For this purpose, the simulation was performed using a 
special procedure, in which in all collisions before and af- 
ter the emission of the shower, the ion deflection was can- 
celed (the single-collision model with impact-parameter- 
dependent energy losses). These results are also shown 
in Fig. [4}d as a dashed curve and demonstrate that the 
multiple and plural scattering events themselves are of 
minor importance for the formation of the surface peak. 
In summary we can conclude that the surface peak re- 
flects mainly the correlation between scattering and en- 
ergy loss: at this low energies ions scattered in deeper 
layers have more chances to exit from the sample if they 
do not re-scatter strongly on atoms of the upper lay- 
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FIG. 5: Results of the simulation of 100 keV He + ions scat- 
tering from a Si crystal. The scattering yield is accumulated 
from depth ranges of (a) 5-30 nm and (b) 30-55 nm. The dis- 
played angular range lies around the <112> crystal direction. 



ers and, consequently, the energy losses in the passage 
through these layers are also small. 

It is worthwhile to note that the program TRIC passes 
here a serious test: the simulation with shower generation 
produces results identical to those obtained when running 
the program in the direct simulation mode. 

Turning now to medium energy ion scattering (MEIS), 
we show in Fig.[5]thc angular distribution of 100 keV Hc + 
ions scattered on a Si(100) crystal. The geometry (the 
<112> axis is in the center of the position-sensitive de- 
tector) and the two depth ranges chosen for collection 
of the scattering events are equivalent to the experimen- 
tal conditions used by Kobayashi [lj5[. The results of 
the simulations are also very similar to the experimental 
results: even at such relatively small depths the washing- 
out of the blocking pattern is well seen, first of all for the 



8 




o Si 
• Y 



[1 00]/[1 1 1 ] 




40 45 50 55 60 65 70 



1.2 



V) 
h — ; 

p 
>- 



[011]/[100] 




0.6 



70 75 80 85 90 95 100 

Scattering angle, dg. 

FIG. 6: Comparison of our simulations (solid curves) of an- 
gular scans for 100 keV p scattering from a monolayer-thick 
YSi2 film on Si(lll) with the result of a VEGAS simulation 
(dashed curves) and experimental data [2(| (dots). The direc- 
tion of the incident beam is parallel to (a) the <I00> and (b) 
<011> directions of the Si substrate while the detector lies in 
the plane (011). The structure of surface layers is shown at 
the top. The yield is normalized to the Rutherford cross sec- 
tion. The dash-dotted curve in the lower graph shows results 
for a trial structure with Y atoms located at the top layer. 



narrow channels. Although this effect of re-channeling 
is well predictable, its detailed demonstration in the re- 
ferred experiment is rather interesting and our simula- 
tions support these results. Notice that, with 10 5 ions 
being sent to the crystal in the simulation, the total yield 
over the detector amounts here to ~ 0.01 only (remind 
that, in the direct simulation, this would be an expected 
number of counts). In this sense, the results of our sim- 
ulation are unique, it is hardly possible to obtain them 
using other known methods. 

In Fig. [6] we show MEIS angular scans for the yield 
of 100 keV protons scattered from Y atoms of a YSi2 
monolayer cpitaxially grown on Si(lll) (a side-view of 
the structure is shown in the picture) . The results of the 



measurements and of the simulations with the program 
VEGAS are taken from Ref. [2(J. The dips in the scans 
are due to blocking of the scattering from Y atoms by 
Si atoms of the upper layers. All parameters in the two 
simulations are taken to be identical. The achievement 
of a perfect agreement by optimization of parameters of 
crystal structure and lattice dynamics is, in fact, the ac- 
tual goal of the MEIS analysis [2(|. Here, we emphasize 
only that results of our simulation coincide well with the 
results of the VEGAS simulation. In principle, this is 
predictable for such thin layers (see Introduction). As 
an additional test of our simulation procedure we per- 
formed simulations also for the case when the top layer 
is terminated by Y atoms instead of Si atoms. As is seen 
from the figure, the scattering yield in this case simply 
reproduces the Rutherford cross section. 

The most difficult problem for simulations is the de- 
scription of conditions where plural scattering is strongly 
effective. Besides low energies, this happens also in the 
medium and high-energy cases where heavy-ion beams 
and/or heavy-atom targets are considered. For the case 
of scattering from amorphous samples, Biersack et al. (2lj 
proposed the fast version of the code TRIM [TtJ where 
the consumption of computer power is reduced at the cost 
of additional approximations. Figure [JJl shows their re- 
sults for scattering of 100 keV protons from a 1000 A- 
thick gold foil, together with the results of simulations 
done with the shower approach. As seen, the two ap- 
proaches produce in fact identical results, also for similar 
computation efforts. The only visible difference is at the 
low-energy tail where yield is entirely due to the plural 
scattering. Therefore, these results can be regarded as a 
confirmation of the accuracy of the TRIM approach by 
comparison with an exact treatment of the model with a 
similar structure of the target. 

In Fig. [TJj the simulation results are compared with 
results of an experiment [23| of 280 keV a-particles scat- 
tering from a 1130 A-thick Pt foil. In general, the agree- 
ment is satisfactory also here. In the two cases shown in 
Fig. [7J the shape of the spectra differs significantly from 
those predicted by the single-collision model: the yield 
is higher for the main part of the spectra, also the low- 
energy tail which is entirely due to the plural scattering 
cannot be predicted at all by the single-collision model. 
To reproduce all the features above the fluctuation level, 
the showers must be generated in wide cones. In the 
considered cases the double-cone approach (see Sect. II) 
was used with a half-width of 40° for the internal cone 
and of 160° for the outer cone. In this way, at least the 
double-scattering events are treated with special efforts. 
The effect of the usage of the two cones is illustrated 
in Fig. [JJi where, for comparison, the simulation results 
obtained in the single-cone approach are included. The 
level of fluctuations in this case suggests that much larger 
efforts are necessary to achieve the same level of preci- 
sion of simulation as in the double-cone approach. To 
demonstrate the role of plural scattering in more detail, 
we show (Fig. [JJd) separately the contribution of showers 
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FIG. 7: Backscattering energy spectra of 0.1 MeV p scattered 
from a 1000 A Au film (a) and of 0.28 MeV He ion scattering 
from a 1130 A Pt foil (b). The simulation results and the cal- 
culations in the single-scattering approximation are compared 
with the RBSTRIM simulation [22] and with experimental 
data [23|]. Results of simulation in the single-cone approxima- 
tion are shown in the top panel. The lower graph shows the 
partial contributions of ions from the internal (MS) and outer 
(PS) cones. 



emitted by ions moving in the outer cone. It is seen that 
these histories of collisions completely explain the nature 
of the low-energy tail and, in general, contain mainly the 
effect of plural scattering. On the other hand, the par- 
tial spectrum of ions in the showers produced directly 
by the primary ion is influenced by multiple scattering. 
It looks surprising at a first glance that, in contrast with 
the predicted increase of the yield compared to the single- 
scattering spectrum [U, [25[ , this partial spectrum shows 
the inverse ordering due to multiple scattering. This is, 
however, natural because in this partial spectrum we do 
not consider the compensation of the transport out of the 
internal cone by the transport coming from the external 
region. 

As a last example we reproduce in our simulations the 
experiment [2rj| performed in the early times of chan- 
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FIG. 8: Channeling and blocking wells calculated for the case 
of irradiation with 1 MeV protons of a tungsten crystal along 
the <100> axial direction with detection at a 135° scattering 
angle and in the inverse geometry. The yield is calculated in 
the energy window corresponding to scattering from a depth 
range between 2500 and 3500 A. The blocking dip normalized 
to the thickness of the layer along the beam direction is also 
shown. 



neling studies with the special aim to test the Lindhard 
time-reversibility rule [10| . Here a proton beam of high 
energy (1 MeV) is incident on a W crystal along the 
<100> direction (close to the surface normal) and the 
protons scattered by 135° into a random direction of the 
crystal lattice are detected. In the inverse situation the 
beam and detection directions were interchanged. The 
angular spread in the beam and the aperture of the de- 
tector were approximately equal, 0.1°, and the yield was 
measured in an energy window corresponding to a layer 
of 1000 A thickness at a depth of 3000 A. The yield was 
measured as a function of the beam misoricntation in 
the first case and as a function of the detection angle 
relative to the <100> channel in the second. The simu- 
lated channeling and blocking wells are shown in Fig. [5] 
Their widths are similar and coincide well with the ex- 
perimental results. To achieve also the coincidence of 
the absolute values we had to normalize the yield to the 
path length of the incoming ions inside the considered 
layer (by multiplying the yield by cos &i n in the second 
case, where @i„ is the angle of beam incidence relative 
to the surface normal). The meaning of this procedure 
is explained in the next section. In general, the time- 
reversibility is confirmed by the simulation similarly as 
in the referred experiment. 



IV. DISCUSSION 

In this section we discuss the advantages of the pro- 
posed approach in comparison with those used currently 
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in simulations of swift ion-solid interaction. First, the 
shower approach is equivalent in accuracy with the direct 
simulation as performed by the code MARLOWE. The 
difference lies only in the much higher efficiency (lower 
level of statistical fluctuations) as illustrated in Sect. II 
and III. This is in fact of primary importance since, in 
many situations, extensive analysis of experimental data 
becomes feasible. Furthermore, for medium and higher 
ion energies a proper simulation of multiple and plural 
scattering is not possible at all by other methods. For 
medium-energy ion scattering from amorphous samples, 
the accelerated version of the code TRIM [2l| seems com- 
petitive but this is achieved at the cost of additional ap- 
proximations. 

Barrett's approach followed by later developments 
like the programs FLUX @ and UPIC Q, is based on 
the single-scattering model and has consequently its re- 
gion of applicability restricted to the cases where mul- 
tiple and plural scattering can be neglected. In princi- 
ple, the calculation of the close-encounter probabilities in 
this approach bears some similarity with our procedure 
of shower generation. In this way, the picture of colli- 
sions with small impact parameters is well reproduced. 
However, only one of the two segments of the trajectory 
is described realistically while the other is approximated 
by a straight line. 

The reversing procedure used in the programs VE- 
GAS @ and MATCH Q3 is closely related with the 
shower approach and it is rather instructive to compare 
both strategies in more detail. First, to fill some gaps 
in the published descriptions, we present here a general 
analysis of the basis of this approach in a systematic way. 
The problem is formulated in terms of the fluxes of in- 
coming and outgoing ions properly convoluted at certain 
lattice sites. Analogously as for the incoming flux, a cer- 
tain distribution of initial conditions (see below) has to 
be assumed for the calculation of the time-reversed out- 
going flux. Therefore, the problem is formulated from 
the outset differently than in the direct simulation: the 
simulation is aimed to obtain the integral yield in some 
range of exit directions Afl and of energies of scattered 
ions E -f- E + AE (the detector acceptance) . In order to 
treat the problem in the most general way, we will not 
assume in the following that Ail or AE are small. 

Further, as in any non-direct approach, special care 
must be taken in order to ensure that the statistics of 
thermal displacements of target atoms is treated cor- 
rectly in order to guarantee that, in this respect, the ap- 
proach is equivalent to the direct simulation. Basically, 
since the ion velocity is much larger than the thermal 
velocities of the target atoms, each ion passing through 
the sample volume sees a frozen pattern of atom dis- 
placements. The configurations of displacements must 
be chosen randomly according to the model assumed for 
the description of thermal vibrations. In the reversing 
approach, as also in the shower approach, we use the fol- 
lowing circumstance (27j : when averaging the scattering 
yield over different configurations of atom displacements 



one may replace the contribution of each configuration by 
the average over the displacements of one of the atoms, 
arbitrarily chosen. Note that the probabilities of dis- 
placements of this atom can be dependent on the actual 
(fixed) positions of other atoms; in fact, accounting for 
correlations in the thermal vibrations is needed in any 
realistic description. The possibility of averaging can be 
used to smooth the statistics of scattering events in the 
simulation: the smoothing is achieved by multiple sam- 
pling of the displacements of an atom met along the ion 
path. For every resulting trajectory the same procedure 
can be iterated at the next lattice site (the atom positions 
chosen in the previous collisions must be considered as a 
modification of the initial configuration of the surround- 
ing). Note that the direct simulation is usually performed 
in the same manner, the only special feature here is that 
only one sample is taken at each lattice site. 

Multiple sampling of atom positions allows to obtain 
several scattering events for every impinging ion. The 
advantage of this is that a significant part of the calcula- 
tions is common to all trajectories. A dramatic increase 
of efficiency in simulation using the shower approach is 
achieved due to a special way of sampling which enhances 
the probability to finally reach the detector. To apply the 
strategy of importance sampling we appeal, in fact, to the 
single-scattering model assuming that it can be used at 
least as a first approximation. The reversing approach 
is based on a different idea. The random sampling of 
atom displacements has, in fact, the only purpose to find 
those positions of the atom which result in useful scat- 
tering events. Every such position defines simultaneously 
the corresponding outgoing trajectory which starts from 
the vicinity of the considered lattice site and ends up 
in the detector. Therefore, the problem could be solved 
in the reversed order, starting with the identification of 
these trajectories. A set of such trajectories can be gen- 
erated considering the time-reversed motion of ions start- 
ing from the detector. With this information the random 
sampling of atom displacements is not needed, the rel- 
evant positions of the atom are those which provide a 
connection of the trajectory of the ingoing ion with one 
of the outgoing trajectories. The trajectories are to be 
connected as a result of the scattering on this atom. 

To recognize clearly the relation between the shower 
and reversing approaches we can formally write the con- 
sidered contribution to the scattering yield Yi (i is the 
site number) as 

Yi =<d\T 2 S t T 1 \b>, (2) 

where \b > denotes the flux of ions in the beam before 
entering into the sample volume, Si the operator of scat- 
tering on the atom and operators T\ and T 2 describe the 
transformations of the flux before and after the collision 
with the «-th atom, respectively. In the shower approach 
the resulting transformation of \b > is projected on the 
detector acceptance profile \d > (\b > and \d > can be 
considered as elements of a real Hilbert state). The re- 
versing approach consists in this picture in the action of 
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the two operators T\ and Si followed by the projection 
on the state T^\d >• 

To specify the way this plan can be implemented, we 
need to represent the action of the operator Si explicitly. 
Denoting by uj\ the phase variables of the incoming ion 
entering the lattice site volume, we can write the prob- 
ability that further motion will end up in the detector 



as 



fi(wi) 



dP sc 

dbJ2— Pout- 



(3) 



The integration is performed over the phase variables of 
the ion after the collision L02, {dP sc / duj2)duj2 is the prob- 
ability of scattering into one of the states within du>2 ; this 
probability is related to the probability that the atom is 
located at the appropriate position. Finally, P ou t{^2) is 
the probability that, in the course of its further motion, 
the ion in a state within di02 will leave the crystal with an 
energy and in a direction within the detector acceptance. 

The probability dP sc /du>2 is explicitly determined in 
the case that both the scattering angle 0(b) and the 
energy loss AE(h) in the binary ion-atom collision are 
uniquely determined by the impact parameter b. To de- 
rive the corresponding expression we describe the states 
lo\ and CO2 in a local coordinate system where the z axis 
is aligned with the ion velocity before the collision. As a 
result, 



dP st 

du>2 



<% 2 - yi)S(E 2 - E 1 + AE(b)) 



d 3 R 



dx2d 2 n 2 



£>(R), 
(4) 

where x and y are the coordinates in the scattering plane 
and in the perpendicular direction, respectively. The first 
J-function satisfies the condition that the two trajectories 
must intersect and the second takes into account the re- 
lation between the energies before and after the collision. 
In the case of potential scattering the state of ion motion 
after the collision (given by the coordinate in the scat- 
tering plane X2 and the velocity direction 112) is uniquely 
determined by the atom position R. The Jacobian of this 
relation appearing in J3]) is expressed as 



d 3 R 



dx2d 2 ri2 



sin 2 9 



db 
dB 



1 da 
sin 9 dfl ' 



(5) 



where b is the impact parameter corresponding to the 
scattering angle 9 and da / d£l is the scattering cross sec- 
tion. Finally, D(H) in (01 is the density of distribution 
of thermal displacements of the target atom. Note that 
R, the atom position which results in the considered col- 
lision, is a function of u>\ and u)2- 

To confirm the validity of Eqs.(j4]) and ([5]), let us calcu- 
late the angular dependence of the probability of scatter- 
ing on a single atom (the variation due to the uncertainty 
of the atom position R). First, we integrate both sides 
of Eq.(U) over y 2 and E 2 and this cancels the two 6- 
functions. The coordinate X2 is related to the coordinate 
z c of the point of crossing of the scattering asymptotes, 



in such a way that dx2 = sin9<i2 c . In addition, we take 
into account that, for a given scattering angle, z c is di- 
rectly related to the coordinate z of the atom: dz c = dz. 
As a result we arrive at the familiar expression: 



d 3 P x , 



da 



dzd 2 vi2 dfl 



D(R) 



(6) 



which expresses nothing but the concept of differential 
cross section da/dVL. 

The contribution of scattering from the considered 
atom to the yield at the detector is obtained as an inte- 
gral of the probability © multiplied by the flux 
of ions of the beam reaching this lattice site: 



Y, 



dwi$> b Pi. 



(7) 



Further, the variables u>\ and u>2 are uniquely related 
to uji n and oj ou t respectively, the parameters of motion 
of the ion when it enters (exits) the sample volume. It 
is useful to refer directly to the latter parameters and, 
for this purpose, we replace the integrations over oj\ and 
(x>2 in ([3]) and (J7J) by integrations over uji n and u! ou t, re- 
spectively. Then, the Jacobians J\ = \duj\/ duji n \ and 
J2 = \diL>2/dio ou t\ must be added into the integrands. 
Note only that the relation uj\{oJi n ) cannot be defined for 
all ujin (some initial conditions uji n result in trajectories 
which never pass through the vicinity of the considered 
lattice site). The same is possible for the pair of vari- 
ables UJ2 and ui ou t- To account for this, we assume the 
value of the respective Jacobian in such cases to be zero. 
Now, combining the above equations, we arrive at the 
expression: 

Yi = J duii n Ji n <&b J duj out JoutPouf 

■ S(y 2 - V1ME2 -E 1+ AE(b))^—^-D{R). (8) 

sin 9 di I 

As a function of io ou t the probability P ou t is easily evalu- 
ated: P ou t = 1 if uJout is within the region of acceptance 
of the detector, otherwise P ou t = 0. One can interpret 
this as the operation of projection of the outgoing flux 
on the exit states, the set determined by the detector 
acceptance; the "density of states" is here uniform. In 
general, all aspects of Eq. (J2J) arc reproduced in the above 
equation. Particularly, the operators 7[ and 7~2 are rep- 
resented as the transformations of variables u>i n (u>i) and 
(•<-W (o;2 ) with the corresponding Jacobians. 

Eqs. (J2J) and (JSj> show that the two segments of the 
ion trajectory, before and after the collision with the 
considered atom, can be treated analogously. The flux 
7^ _1 |d > can be represented as the result of the time- 
reversed evolution of the "beam" with a uniform profile 
(uniform distribution of angles and energies within the 
detector acceptance and, also, uniform distribution of 
initial coordinates at the sample surface). This shows 
a way in which the integral © can be evaluated using 
the Monte-Carlo approach. The relations of variables 
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wi(wi n ) and uj^^out) can be reconstructed by the calcu- 
lation of a number of trajectories of incoming ions with 
the corresponding distribution of initial conditions and 
a number of trajectories of outgoing ions calculated in 
time-reversed mode (including negative energy losses). 
The densities of the two fluxes approaching the consid- 
ered lattice site are represented by the factors Jm^b and 
JoutPout-, respectively, in the integrand of |5]), 

The integration itself, however, presents here a serious 
problem. First, it is difficult in practice to reconstruct 
in numerical form the distribution of the fluxes in the 5- 
dimensional section of the phase space. And, second, the 
presence of ^-functions in ([8]) results in the need of some 
special procedure. For this reason, the standard numeri- 
cal integration is practicable only in a simplified version 
of this approach, as performed by the computer code VE- 
GAS, where the energy is not resolved and the scattering 
angle is assumed to be fixed. In this case the procedure is 
reduced to the convolution of two two-dimensional fluxes 
at a closed surface containing the considered atom. 

If the fluxes are represented by sets of trajectories, 
the conditions expressed in jS]) by the (^-functions (the 
"conditions of crossing" for a pair of ingoing and outgo- 
ing trajectories) cannot be satisfied exactly. Therefore, 
to accumulate a representative statistics, one has to in- 
troduce in these conditions some tolerances, sacrificing 
thereby the accuracy of the simulation. In the numerical 
integration where the density in phase space is presented 
by histograms, the sizes of the bins of the histograms 
play the role of the tolerances. Another procedure for 
the integration is implemented in the program MATCH. 
Pairs of incoming and outgoing trajectories are matched 
to identify that crossing at the considered lattice site. 
The tolerances for coincidence of coordinate y and for 
the energy relation, t v and te, are specified beforehand 
and, in terms of Eq.®, the whole procedure can be for- 
mulated as follows. The tolerances are introduced as a 
replacement of the ^-functions by normalized pulse func- 
tions 



6{t) => U(t) = 



1 



1 if \t/r\ < 1/2 
if \t/r\ > 1/2 



(9) 



The integration over u) ou t in |8]) is associated with the 
Monte-Carlo sum: 

SoMlAE^ 

OJjJ ou tJout^out =? 77 2_j ' V / 

1* out . 

%—\ 

where So is the area of the unit cell of the crystal surface 
and N out is the number of calculated outgoing trajecto- 
ries. To obtain the yield of scattering for Ni n ions in the 
beam, we replace the integral over u>i„ (together with the 
Jacobian Jj„ and the beam profile $(,) simply by a sum 
over all impinging ions. The terms of the resulting dou- 
ble sum are determined by the rest of the integrand in 
©: 



SqAHAE 1 da 
N utTyT E sin 9 dCl 



D(R). 



(11) 



Hence, each crossing of trajectories contributes to the 
accumulated statistics with its own weight. 

The weight (fTl) differs from the weight actually ap- 
plied in the program MATCH, where it is simply taken 
as the product of the cross section and the atomic den- 
sity. The dependence on the angle 6 in (jTTJ) includes, 
however, an additional factor l/sinO. The necessity of 
such factor is immediately recognized when one tries to 
calculate the yield of scattering from one atom ([6]) by the 
simulation. The factor 1 / sin in (|11|) is compensated in 
this case by the number of crossings of trajectories, which 
is proportional to sinO; the latter follows from a simple 
geometrical consideration. The first fractional factor in 
(jlip plays the same role, the normalization on the num- 
ber of crossings. As a result, the Monte Carlo sum, the 
sum of weights (|TTj) , gives the absolute value of the scat- 
tering yield Y^. Interestingly, the first fractional factor 
in the weight pip is reminiscent of the quantum defini- 
tion of the number of states, the volume of phase space 
divided by the Planck's constant. Here, analogously, the 
tolerances r y and te represent the volume per state for 
states of ions exiting from the sample volume; the phase 
space attains thereby a coarse-grained structure. 

As stated above, the described procedure is equivalent 
to the simulation in the direct mode, the result represents 
the contribution of the chosen configuration of atoms dis- 
placements additionally averaged over displacements of 
the i-th atom. As stated above, such averaging does not 
bias the results of the simulation and it might be ex- 
pected that the resulting smoothing of the accumulated 
statistics will improve significantly the efficiency of the 
simulation. However, one can meet serious problems in 
the practical implementation of this approach. The first 
problem is the choice of reasonable values for the tol- 
erances T y and te- In fact, the use of tolerances intro- 
duces some inaccuracy into the description of scattering 
and, consequently, the tolerances must be small to pro- 
vide sufficient accuracy of the simulation results. On the 
other hand, small tolerances imply smaller probabilities 
of trajectory crossings and, therefore, the efficiency of the 
calculation is lowered. Thus, the optimal values are the 
maximal values of t v and te which still provide sufficient 
accuracy. 

The inaccuracies have the following origin. Replace- 
ment of the 5-functions in j8]) by finite-width pulse func- 
tions is allowed only when the fluxes Ji n and J ou t are 
sufficiently smooth. In fact, we assume here that each cal- 
culated trajectory of ion motion can be associated with a 
family of close trajectories, the family limited by the tol- 
erances, and they are treated simultaneously. The criti- 
cal points in this assumption are the following. First, the 
simultaneous treatment of the family of trajectories as- 
sumes that they all reach the detector. This means that 
the tolerances should correlate with the detector accep- 
tance, the "resolution" in the connection of trajectories, 
so to say, must be finer than that of the detector. In this 
connection, one should also take into account that, af- 
ter many collisions with atoms on the outgoing path, the 



13 



ions of the family can exit from the sample in strongly 
different directions and, as a result, not all of them will 
be detected. It is clear that this difficulty increases when 
the ion energy is increased (scattering from large depths) . 

Applying this approach one must decide also at which 
lattice site the desired matching procedure will be per- 
formed. Assuming, as in Eq.©, that the site is the same 
for all incoming ions, one has to choose the site in the 
top layer, the first one met along the ion path. Other- 
wise, this atom should be considered as being at a fixed 
position and, therefore, direct scattering in the direction 
to the detector would be possible. However, in such im- 
plementation of the approach, the problem of strong sta- 
tistical fluctuations is not entirely eliminated. The fluc- 
tuations manifest themselves as extremely large values of 
the weight (flT) when two trajectories intersect at a small 
angle O. Such fluctuations rarely appear in the cases 
when a large-angle scattering (with small cross section) 
happens at some other lattice site. Here, such scattering 
events (happening at larger depths) are present in the 
outgoing trajectories. In fact, this means that the large- 
angle scattering is treated here in a way equivalent as in 
the direct simulation and this is a clear disadvantage. 

The procedure can be modified in such a way that the 
connection of ingoing and outgoing trajectories is per- 
formed only at the points of large-angle scattering. The 
modification consists in the following. First, the inter- 
esting site will be chosen as specific for each ingoing ion. 
The whole trajectory of ion motion in the frozen lattice 
is inspected to find the site where the ion most closely 
approaches the equilibrium position of the atom. If the 
matching procedure is performed at this site, the range 
of impact parameters of possible collisions is restricted, 
it is limited by the amplitude of thermal vibrations of 
the atom. In the case of low-energy heavy ions (strong 
interaction) this means an effective restriction of the pos- 
sibility of small-angle scattering. However, even at low 
energies, this does not ensure an appropriate low level of 
fluctuations of the weights Wi. At higher energies this be- 
comes impossible at all, the allowed values of the impact 
parameter correspond to large variations of 0. 

It is readily recognized that this is the same problem 
as that solved by the shower approach. Again, the main 
difficulties lie in the treatment of the plural scattering. 
To solve this problem we can use again the main idea of 
the shower approach, which is to emphasize the events 
of scattering in the direction to the detector. The only 
difference will be that the multiple sampling of atom dis- 
placements from the "hot" region is replaced by a corre- 
sponding selection of crossings of trajectories. Namely, 
the position of the atom necessary for the connection of 
ingoing and outgoing trajectories must be located within 
the "hot" region. 

Therefore, the reversing approach turns out to be noth- 
ing but a special version of the shower approach. In or- 
der to choose, for a given application, which of the two 
versions is more convenient, one should consider the fol- 
lowing. First, the simulation with the reversing version 



cannot be considered as an exact treatment inside the 
binary collision model since the finite tolerances for con- 
nection of trajectories represent an unavoidable source 
of inaccuracies. Practically, the only way to control the 
inaccuracies is the repetition of the simulation using dif- 
ferent values of t v and rg. Obviously, this is a serious 
drawback of the approach. Another drawback is that the 
whole procedure of the reversing approach is too cum- 
bersome for an extensive usage and for developments of 
the model of interaction. The arguments above show 
that, with some simplifications, this approach can be ap- 
plied to the simulation of scattering of low-energy heavy 
ions (as it was used before) but is hardly practicable at 
higher energies. It is even difficult to estimate, for exam- 
ple, what kind of efforts would be necessary to reproduce 
the data presented in Fig. [7] For the program TRIC, this 
is also a difficult problem, the calculations with an ordi- 
nary PC take about 20 hours. One can conclude, how- 
ever, that in general, this program is capable to address 
by simulations a wide class of problems of scattering and 
recoil emission. 

To finish this section we make some remarks intended 
to clarify possible consequences of the reversibility rule 
for the interpretation of experimental results. It is seen 
in Eq.© that, under the assumption of a purely poten- 
tial motion of the ions (the Jacobians Ji n = J ou t = 1), 
the yield of scattering from one lattice site is symmetric 
under interchange of the beam and detector directions 
provided that the flux in the beam is also uniform and 
the respective phase-space volumes are equal. The lat- 
ter condition is less relevant because the difference can 
be simply accounted for by a factor in the yield. There- 
fore, we can conclude that the yield of scattering from 
one atom for a given beam-detector configuration is de- 
termined if it is known for the inverse situation. How- 
ever, in the measurements of the yield in a certain energy 
window, as ordinarily made, the effective number of con- 
tributing atoms can be different. This fact was taken 
into account in the transformation of the data shown in 
Fig.! 



V. CONCLUSION 

The shower approach proposed in this paper solves ef- 
fectively the main problem of simulation of ion scatter- 
ing from solids within the binary collision model, which 
is elimination of the violent statistical fluctuations in the 
Monte-Carlo sum. This is achieved as a result of specific 
effective improvements of the direct simulation approach. 
The improvements consist in the use of the strategies of 
importance- and stratified sampling, special ingredients 
of the Monte-Carlo method. As a result, the computer 
power required for simulation is reduced by several or- 
ders of magnitude. This is, in fact, a decisive advantage 
allowing to address simulation problems which cannot be 
treated with other methods. As discussed in Sect. IV, our 
method avoids also many shortcomings inherent to alter- 
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native approaches. It is argued in particular that, in fact 
for the first time, this method allows a reliable treatment 
of the rare events of plural scattering. Such possibility 
is especially important because the plural scattering is 
also not amenable to theoretical treatment and the only 
proposed method of simulation, the modified version of 
the code TRIM, has significant restrictions. First of all, 
it is capable to treat only amorphous structures. 

We performed a detailed analysis of the basis of the 
reversing approach, as applied in the programs VEGAS 
and MATCH. The latter program offers the only alterna- 
tive for an exact treatment of the binary collision model, 
including multiple and plural scattering. It is argued in 
particular that the corresponding procedure for match- 
ing must use a similar scheme of importance sampling as 
used in the shower approach. The main difference lies in 
the method of sampling. Basically, if the scattering from 
shallow depths is considered, the two methods are equiva- 
lent. The applicability of the reversing approach at larger 
depths needs, however, a special justification. Also, right 
from the beginning, this approach uses an approximation 
(the finite tolerances assumed for the connection of the 
two segments of the ion trajectory). Note also that the 
procedure used in its practical application is much more 
cumbersome than that used in the shower approach. In 
fact, in simulations for medium and high energy ions this 
approach becomes impracticable at all. 

Simulations with the code TRIC can be performed for 
large classes of crystal structures of the samples and pro- 
vide a detailed picture of scattering or recoil yields in 
the form of energy and angular distributions. All these 
qualities are promising for a wide use of the developed 



computer code both in basic research and in the analysis 
of materials. In particular, this provides the possibility to 
efficiently compare measurements with simulations made 
for many trial structures allowing in this way a precise 
structural analysis. 

This paper shows several examples of the use of the 
program although, of course, it is not possible to cover all 
potential applications (e.g. simulations of the sputtering 
or total reflection yields, well possible in this approach, 
are also interesting applications of the code TRIC). The 
illustration examples in Sect. Ill are chosen to demon- 
strate the capabilities to solve specific problems and to 
show the accuracy of this method in comparison with 
others. In particular, we show the capability to simu- 
late the interaction of different ions of low and medium 
energies with solid matter including complex structures, 
to calculate the yield of scattered ions and recoils and 
to reproduce their energy spectra and angular distribu- 
tions. In addition, we address the interpretation of the 
time reversibility rule and provide additional insight into 
the origin of the surface peak. 
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